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The epidermal growth factor receptor (EGFR) signaling network is activated in most solid tumors, 
and small-molecule drugs targeting this network are increasingly available. However, often only 
specific combinations of inhibitors are effective. Therefore, the prediction of potent combinatorial 
treatments is a major challenge in targeted cancer therapy. In this study, we demonstrate how a 
model-based evaluation of signaling data can assist in finding the most suitable treatment 
combination. We generated a perturbation data set by monitoring the response of RAS/PI3K 
signaling to combined stimulations and inhibitions in a panel of colorectal cancer cell lines, which 
we analyzed using mathematical models. We detected that a negative feedback involving EGFR 
mediates strong cross talk from ERK to AKT. Consequently, when inhibiting MAPK, AKT activity is 
increased in an EGFR-dependent manner. Using the model, we predict that in contrast to single 
inhibition, combined inactivation of MEK and EGFR could inactivate both endpoints of RAS, ERK 
and AKT. We further could demonstrate that this combination blocked cell growth in BRAF- as well 
as KRAS-mutated tumor cells, which we confirmed using a xenograft model. 
Molecular Systems Biology 9: 673; published online 11 June 2013; doi:10.1038/msb.2013.29 
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Introduction 

The signal transduction network downstream of the epidermal 
growth factor receptor (EGFR) has received much attention, as 
a majority of human cancers shows mutations leading to 
hyperactivation of the network (Hanahan and Weinberg, 
2011). Based on detailed mechanistic understanding of the 
network, a large number of targeted therapies has been 
developed (Herbst et aZ, 2004; Roberts and Der, 2007; Prenen 
et aU 2010). However, despite positive treatment responses in 
some patients, a large fraction of patients do not benefit even if 
molecular markers such as KRAS or BRAF mutation status are 
used to stratify patient groups (Karapetis et aU 2008; Walther 
et aU 2009; Roth et aU 2010). 

One reason for the somewhat disappointing response rate to 
these therapies is that they have been developed using the 
concept of linear signaling pathways downstream of the 
receptor. However, the EGFR signal is propagated through a 
complex network (Bublil and Yarden, 2007), involving cross 
talks to parallel pathways (Porter and Vaillancourt, 1998} and 
strong feedback loops on different levels (Bliithgen and 



Legewie, 2008; Legewie et al 2008; Cirit et al 2010; 
Avraham and Yarden, 2011). Quantitative analysis of these 
regulatory principles suggested that strong feedbacks can 
neutralize drug treatment (Friday et aU 2008; Cirit et aU 2010; 
Sturm et al 2010; Fritsche-Guenther et al 2011). 

Mathematical modeling of signaling networks can help to 
understand the behavior of these complex networks, and can 
be used to simulate the effect of drugs in such a network. The 
structure of these mathematical models can be directly 
deduced from pathway maps (Oda et al 2005). Detailed 
mechanistic models based on Ordinary differential equations 
(ODE) have been developed for the EGFR signaling network 
(Kholodenko et al 1999; Schoeberl et al 2002; Nelander et al 
2008). However, for such detailed models the parameteriza- 
tion remains a major challenge. More coarse-grain modeling 
approaches, such as logical models or non-mechanistic 
statistical models require less data for parameterization 
(Kreeger et al 2009; Morris et al 2011; Saez-Rodriguez et al 
2011, 2009; Tentner et al 2012). These approaches allow 
qualitative predictions, but typically fail to deal with feedback 
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loops or do not provide mechanistic insights. The approach we 
chose for this study is termed modular response analysis 
(MRA), which resides between the qualitative nature of 
Boolean models and detailed mechanistic models. It provides 
a framework to calculate the response of a linear approxima- 
tion of an ordinary differential equation model to a perturba- 
tion (Bruggeman et aU 2002; Kholodenko et aU 2002), and has 
been developed to discover and parameterize networks from 
systematic perturbation studies (Santos et aU 2007; Stelniec- 
Klotz et aU 2012). The parameters of an MRA model are so- 
called local response coefficients that quantify how strong a 
change in activity of one node directly affects the activity of 
another node. These models then allow to quantitatively 
analyze feedback regulation, feedforward loops as well as 
cross talks, which is of major interest as these network motifs 
have major effects on drug sensitivity and network behavior 
(Friday et aU 2008; Cirit et al 2010; Sturm et al 2010; Fritsche- 
Guenther et al 2011). 

In this work, we exposed a panel of colon cancer cell lines to 
different stimuli and pharmaceutical inhibitors, and measured 
key signaling molecules in a medium-throughput approach. 
The data generated by this approach were then used to 
parameterize MRA-based mathematical models, which gener- 
ated quantitative maps of the wiring between signaling 
molecules. We focused our efforts on RAS-mediated signal 
transduction pathways, as they are currently in the strategic 
focus of targeted therapeutics in solid cancers. We were able to 
identify feedbacks and cross talks of therapeutic relevance. 
Our model predicted that EGFR-directed therapeutics might be 
effective even in tumors carrying a mutation in RAS, if they are 
provided in combination with RAF or MEK inhibitors. We 
confirmed our predictions by phenotypic assays and a 
xenograft model. 
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Results 

A pipeline to model signal transduction networks 
in cancer cell line panels 

We developed a combined experimental and theoretical 
approach to dissect signaling networks in cancer cell lines 
to generate predictive mathematical models for their signal 
transduction pathways. The workflow of our pipeline 
was as follows (Figure 1): we selected a panel of six colon 
cancer cells that were genotyped for common oncogenes using 
targeted next-generation sequencing. Subsequently, the sig- 
naling network in these cell lines was perturbed using small 
molecule inhibitors and growth factors, and combinations 
thereof. Before and after perturbation, phosphorylation of key 
signaling molecules was quantified. This high-dimensional 
data set was then used to parameterize mathematical models 
of the signaling events. These models were simulated to 
predict effects of inhibition and potential combinatorial 
treatment. 



Network quantification by a systematic 
perturbation screen 

We first chose six colon cancer cell lines for our systematic 
perturbation screen, LIM1215, HCT116, SW403, SW480, HT29 



Simulation and 
experimental verification 



Figure 1 General outline of the study. A panel of six colon cancer cell lines was 
chosen, and profiled by sequencing selected cancer-related genes. The cells 
were then systematically perturbed with four kinase inhibitors and two ligands of 
growth receptors, and phosphorylation of key signaling proteins was measured 
using the Luminex platform. These data were used for parameterizing a 
mathematical model, which was then used to predict combinatorial treatments. 



and RKO. We used targeted sequencing of 46 genes to verify 
that these cells represent a panel that reflects the genetic 
diversity of colon cancer (see Table 1 and Supplementary 
Table SI). 

To generate a semi-quantitative database for further 
mathematical modeling, we decided to measure phosphoryla- 
tion of selected signaling molecules for combinatorial pertur- 
bations. In line with previous approaches (Nelander et al, 
2008; Saez-Rodriguez et al 2009, 2011; Morris et al 2011), we 
used a pair-wise design of perturbations, where each inhibi- 
tion is combined with each stimulus. Specifically, we 
stimulated the cells with two growth factors, TGFa and IGF, 
activating the EOF receptor and the IGF-receptor, respectively 
(shown in red in Figure 2 A) . The cells have been pre-incubated 
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Table I Mutation spectrum of the cell line panel 


Gene symbol 


LIM1215 


HCT116 


SW403 


SW480 


HT29 


RKO 


TCGA patients 


ABLl 


P309A 


Y257C 












APC 




K1462R 




A1457T/K1462R 








BRAF 










V600E* 


V600E* 


13 


CTNNBl 


T41A^ 














FGFR3 




S400R 


G12V*'" 


S400R 








KRAS 


A146T 


G13D* 


G12V*'" 






6,8"28 


PIK3CA 




H1047R* 








H1047R* 


4 


SMAD4 










Q311X*'" 






SMO 




V404M 












STKll 




G58S 




G58S 








TP53 






R273H" 


R273H" 


R273H*'" 




26 



Shown are non-silent mutation results from sequencing of known mutated regions in 46 cancer-related genes. SNPs reported by Cancer Genome Project are marked 
with asterisk (*) , and " indicates homozygous mutations. The number of patients found in The Cancer Genome Atlas to harbor this particular mutation (in the order of 
appearance) is shown under TCGA patients. 
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Figure 2 Generation of systematic perturbation data. (A) Perturbations consisted of two ligands (red nodes), and four pharmacological inhibitors (yellow flashes). 
These were applied alone and in inhibitor-ligand combinations for the indicated time points. Then eight phosphorylation signals were measured (blue nodes). (B) Log2 
fold change of phosphorylation in response to the perturbations for the six indicated cell lines. Displayed response range was limited to ± 3.5 (approx. 10-fold). Source 
data for this figure is available on the online supplementary information page. 
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for 1 h with pharmacological inhibitors against the kinases 
MEK, PI3K, IKK or GSK-3a/b (flashes in Figure 2A}. As 
signaling typically displays a strong transient response 
followed by a long-term plateau, we performed time-series 
experiments to determine optimal time points for the experi- 
ments (see Supplementary Figure SI). We find that peak 
transients are within the first 10 min after stimulation, and the 
response to TGFa as well as IGF has reached a plateau at 30 min 
after stimulation. Thus, we chose the 30-min time point for 
further experiments, as the interpretation of MRA requires the 
signaling network to be approximately in steady state. We then 
used the Luminex proteomics platform to measure the 
phosphorylation of eight key signaling proteins (AKT^"^""^, 
ERK2™'^''''\ MEj^^S2i7/s22i^ p70S6K^"2i/S424^ IGF-IR^i^^\ 

GSK-3a/b^^^/^^ IkB-a^^^/^^^ andIRS-1^^^^/^^^^) that are within 
or in close proximity to the stimulated pathways and the 
inhibited kinases (shown in blue in Figure 2A}. 

The resulting data sets are shown in Figure 2B as log2 fold 
changes compared with unperturbed controls. When compar- 
ing the six different cell lines, we found strong response 
patterns to the perturbation in all cells, with the exception of 
RKO cells that show hardly any response and which we 
therefore excluded from further analysis. 

All other cell lines, while clearly responding to the inhibitors 
and stimuh, showed subtle differences in their response 
patterns. For example, IRS-1 phosphorylation changed 
strongly in HCT116 cells at certain treatments, but was not 
altered in SW403 cells. To further infer qualitative and 
quantitative differences in the interactions between the 
signaling nodes and to predict effects of combinatorial 
treatments, we decided to analyze the data using a mathema- 
tical model. 



A model-based approach to analyze the 
perturbation data set 

We developed an algorithm that can unambiguously deter- 
mine network structure and parameters from the data set 
based on MRA. The algorithm extends our previously 
developed methodology (Stelniec-Klotz et aZ, 2012} by 
considering multiple perturbations and unobserved nodes. 
First, the model estimates the response coefficients and 
inhibitor strengths for a literature-based starting network. 
Then edges are iteratively removed if they are not supported by 
the data (tested by a likelihood ratio test), until no further edge 
can be removed (Figure 3 A) . 

One of the major obstacles in comparing the networks 
quantitatively between cells was to find a parameterization of 
the model that can be uniquely determined from the data. For 
example, the individual response coefficients along the path 
from EGFR to the activation of AKT cannot be determined 
individually, as the only node measured is the phosphoryla- 
tion of AKT (Figure 2A}. Thus, we had to re-parameterize the 
MRA model. Elegant rule-based re-parameterization algo- 
rithms (Saez-Rodriguez et aU 2009} could not be apphed, as 
they fail for more complex scenarios, where paths branch, 
feedbacks exist and inhibitors are placed within the path. We 
therefore developed an algorithmic approach that replaces 
unidentifiable parameters by an identifiable combination of 
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parameters. These new parameters can then be reliably 
estimated from the data set and compared between different 
cell lines (see Materials and methods. Supplementary Method 1, 
and Supplementary Figure S2} . 



Model fit requires alteration of the network model 

When we apphed our modeling procedure to the perturbation 
data sets of the five selected cell lines, we obtained fits with 
Chi-squared values between 116 and 470 (Figure 3B, dark 
bars). Our algorithm then pruned the network and could 
remove several links in each cell line without significant 
increase in the Chi-squared values (Figure 3B, blue bars), 
indicating that these links are not important in mediating the 
perturbed signals (see an example for the stepwise reduction 
in Supplementary Figure S2) . 

When we compared the model predictions with the data, we 
found that our literature model failed to explain an increase of 
MEK and ERK phosphorylation after IGF treatment that was 
visible in multiple cells. When checking back with literature, 
we found that IGF-1 can lead to activated RAS by triggering the 
activation chain via IGF-IR, IRS-1 and Grb-2 (Florini et aZ, 
1996). Hence, we included a link from IGF to RAF, which 
improved model performance for all five cell lines significantly 
(red bars in Figure 3 A) . 

Additionally, our model could not account for an increase of 
phospho-ERK level upon treatment with IKK inhibitor 
BMS345541, which was particularly pronounced in HCT116 
and SW403. Interestingly, this increase in ERK phosphoryla- 
tion did not coincide with an increase in MEK levels, thus the 
effect of IKK inhibition on ERK cannot be explained via 
upstream components in the pathway. We thus decided to 
include this potential new interaction in the model and to 
repeat the model selection procedure. The inclusion of a link 
from IKK to ERK improved the fit of the model considerably for 
HCT116 and SW403 cells (see decrease of Chi-squared value in 
Figure 3B orange bars) and showed no improvement for the 
other cell lines. 

The observed increase in ERK phosphorylation after 
treatment with IKK inhibitors may be either due to unspecific 
effects of the chosen IKK inhibitor or due to direct or indirect 
regulation of ERK by IKK or its downstream kinases. We 
therefore decided to characterize this interaction further. We 
found that treatment with the IKK inhibitor BMS345541 
resulted in an increase of phosphorylation of ERK within 1 h 
(Figure 3C, left panel), but treatment with two other IKK 
inhibitors had no effect on ERK phosphorylation, although 
they also blocked phosphorylation of IkBa with similar 
strength at this time point (Figure 3C, right panel). From that 
we concluded that the interaction is likely to be a hitherto 
unknown side effect of the inhibitor. Interestingly, while 
phosphorylation of ERK and its cytoplasmic target p70S6K are 
increased in response to BMS345541 (Figure 3D), ERK's 
nuclear activity seems to be decreased, as expression of 
typical immediate-early target genes of ERK such as EGRl and 
FOS is strongly reduced (Figure 3E). This suggests that 
although ERK phosphorylation is increased, ERK activity 
seems to be redirected toward cytoplasmic targets, and 
treatment with BMS345541 may result in a repression of 
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Figure 3 EGFR modeling required network alteration that led to the discovery of a new effect of BMS345541 on ERK. (A) Scheme of the modeling pipeline from the 
starting network to its parameterization. Systematic perturbation data and a starting network serve as input. The core-fitting routine consists of three steps, and is 
illustrated by a four-node example network. First, non-identifiable parameters are detected and reparameterized. Second, the identifiable parameter combinations are 
fitted to the experimental data. Third, connections not significantly contributing to the fit are removed and subsequently the remaining parameters are refitted. The 
resulting network contains information about the strength, direction and sign of its connections. (B) Chi-squared values of the models trained on data of five cell lines. 
Dark bars indicate values for the initial model. Blue, red and yellow bars show values after model reduction of literature network, with a link from IGFR to RAF, or when 
additionally including a link from IKK to ERK, respectively. (C, D) Log2 fold change in time-series experiments after treatment with IKK inhibitors or solvent controls to 
investigate the IKK-ERK relationship in HCT1 16 cells. (C) IKK inhibitor BMS345541 treatment results in an increased phospho-ERK level, whereas treatment with IKK 
inhibitors PHA408, PS1145 and solvent control (DMSO) results in no increase. The inhibition of IkBa phosphorylation is comparable for all three inhibitors after 1 h 
(Luminex, n= 1). (D) Phosphorylation of p70S6K, a cytoplasmic target of ERK, increases after 60 min of BMS345541 treatment (western blot, n = 2) when compared 
with DMSO and PBS control. (E) Expression of ERK target genes EGR1 and FOS decreases after BMS345541 treatment (qrt-PCR, n = 2). Source data for this figure is 
available on the online supplementary information page. 



typical genetic programs stimulated by ERK. In line with this 
hypothesis, we found that treatment with moderate concen- 
trations of BMS345541 results in impaired prohferation of 
HCT116 cells (Supplementary Figure S3}. Thus, the modeling 
procedure identified an unknown side effect of BMS345541 on 
ERK, and consequently we excluded the data of BMS345541 
from further analysis. 

After these alterations in the network structure, the resulting 
model could mimic most of the responses in quantitative 
detail. Figure 4 A shows the experimental measurement side- 
by-side with the corresponding model fit. For each signaling 
node, measured phosphorylation is displayed in the upper 
row (indicated by filled triangle) and the model simulation in 
the lower row (open triangle). By running the procedure on 
the data with simulated noise added, we confirm that the 
procedure robustly identifies the parameters (Supplementary 
Figure S4). In addition, we confirmed by 100 runs of 
our algorithm that for each data set the model structure 
and parameterization remained identical, irrespective of initial 
parameterization. 
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Model fit uncovers differences in network 
structure 

The topology of the final models is displayed in Figure 4B. 
Edges that our modeling procedure has removed in at least one 
cell line are depicted as dashed lines, with color-coded circles 
indicating the cell line. One interesting difference in the 
topology of the networks resides in the feedback from ERK to 
RAF, which has been removed in HT29 cells. This is in line 
with previous findings that BRAF^^^^^ mutation disables 
MAPK feedback regulation via RAFs (Friday et al, 2008; 
Sturm et aU 2010; Fritsche-Guenther et aU 2011). Furthermore, 
the topology that connects the output kinases to ERK and AKT 
differs. For example, the phosphorylation site measured for 
IRSl is not connected to ERK and AKT in LIM 121 5, SW480 and 
SW403, and only connected to ERK in HT29 cells. Taken 
together, the model-fitting procedure allowed to identify 
qualitative differences in signal transduction networks in our 
cell line panel, some of which can be related to specific 
mutations. 
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Figure 4 IVIodel fit uncoveres differences in network topology and signaling strength. (A) Heat maps for the five cell lines showing log2 fold changes of phosphorylation 
(filled triangles) and the corresponding model simulation (open triangles). (B) Network structure derived from the model fit. Dashed lines indicate edges that are removed 
in those cell lines marked by filled circles (same colors and layout as in A). (C) Model-deduced parameter values representing signaling strength for identifiable paths, 
inhibitors and feedback/cross talk expressed as log2 response coefficients. Orange and blue boxes indicate feedback from ERK via EGFR to MEK, and cross talk from 
ERK via EGFR to AKT. 



Quantitative differences between signaling 
networks 

In addition to qualitative differences in the underlying network 
topology, we were also interested to study how signaling 
networks in these cells differ in quantitative aspects. Thus, 
we decided to inspect differences in the parameter sets. 
Overall, the models had on an average 15 parameters (between 
14 for SW403 and 17 for HCT116}. As these parameters 
typically correspond to comphcated combinations of 
response coefficients, the most intuitive way to inspect these 
parameter combinations was by recomputing these parameter 
combinations such that they correspond to paths in the 
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network. The values of such parameter combinations are 
visualized in Figure 4C. Interestingly, many paths that 
correspond to the intracellular logic, such as the response 
coefficient from MEK to ERK or from ERK to p70S6K have 
comparable values in all cell line models. This suggests that 
certain quantitative aspects of signaling are comparable in 
these cells, despite their heterogeneous genetic background. 
On the other hand, there are also strong differences between 
the cells, most of which relate to external perturbations. 
For example, TGFa does have a strong effect on MEK 
phosphorylation in HT29 cell, but can only weakly activate 
MEK in HCT116 cells. 
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All cells show negative response coefficients for feedback 
regulation from ERK via EGFR back to MEK (orange box in 
Figure 4C}, with strongest negative values in HT29 cells. 



An EGFR-dependent feedback causes cross talk 
between MAPK and AKT signaling 

An interesting aspect of this ERK-EGFR feedback is that it 
connects ERK signaling with AKT activity in an EGFR- 
dependent manner, as ERK inhibits the EGF receptor that also 
stimulates AKT signaling. Thus, inhibition in the ERK pathway 
can lead to an activation of AKT if hgands for the EGFR 
receptors are present. This effect has been previously 
described as a mechanism of drug resistance in tumor cells 
with BRAF mutation (Prahallad et al, 2012). When we 
inspected the response coefficient of the path from ERK to 
AKT via EGFR, we found high negative numbers for HT29, the 
Bi^AF-mutated cell line (blue box in Figure 4C} . However, the 



model fit unveiled that the feedback is present in all cell lines, 
and causes strong cross talk from ERK to AKT in all cell lines, 
including those harboring different KRAS mutations (blue box 
in Figure 4C, Table I). 

To confirm the EGFR-dependent cross talk between ERK and 
AKT, we performed independent experiments in HT29 and 
HCT116 cells, as examples for cells with BRAF and RAS 
mutations. Figure 5 A shows that AKT phosphorylation is only 
slightly increased when a MEK inhibitor is apphed alone, and 
that AKT can be only weakly stimulated with TGFa. In line with 
the cross talk hypothesis, we found that AKT phosphorylation 
was significantly increased (P<0.05} by a factor of 2-3 when 
cells were pre-treated with the MEK inhibitor, confirming that 
the cross talk operates irrespective of mutations in BRAF or 
RAS. We also observed a similar effect when we stimulate with 
EGF for lOmin, another ligand of the EGFR (Figure 5A). 

We then aimed to investigate whether the negative feedback 
directly acts at the level of the EGF receptor, or more generally 
desensitizes growth factor receptors, for example by acting on 
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Figure 5 ERK-AKT cross talk is mediated by EGFR and independent of RAS or H/AF mutation. (A) Increase of AKT phosphorylation in HCT1 16 (/CH/AS mutation) and 
HT29 (BH/^F mutation) cells incubated with MEK inhibitor AZD6244 (0.1 |iM) or its solvent control (DMSO) for 1 h before application of TGFa or BSA for 30 min. Effect 
can be also produced by 10 min EGF treatment in the presence of AZD6244 (1 |iM). (B) Schematic view of the proposed ERK-AKT cross talk via the EGFR feedback 
(red), predicting no effect of MEK inhibition (flash) on stimulation of other growth receptors. (C) Increase of Akt phosphorylation compared with solvent control (DMSO) 
with combinatorial treatment of AZD6244 (5 |iM) and different growth factors (10 min). (D) Response of phospho-ERK and AKT to 10 min EGF treatment for different 
preincubation times of AZD6244 (1 |iM) in HCT1 16 cells. BSA and DMSO are solvent controls for EGF and AZD6244, respectively. All data measured using Luminex 
assays and shown as fold change; significant deviations are indicated. Brackets indicate significant one-sided Mest, P<0.05. Error bars indicate s.d. of 3 samples. 
Source data for this figure is available on the online supplementary information page. 
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shared adaptor molecules (Dhillon et aU 2007) . To disentangle 
the network that mediates the feedback, we stimulated cells 
with the growth factors HGF, FGF and IGF, which do not signal 
via EGFR (Figure 5B). In contrast to EGF, pre-inhibiting either 
cell line with AZD6244 had no significant effect on the response 
to HGF. Similarly, it had no effect on the response to FGF in 
HT29 cells, and only a weak effect in HCT116 cells. When using 
IGF as a control that stimulates primarily ART, pre-treatment 
with MER inhibitors had no effect on ART phosphorylation. 
This suggests that the feedback common to HCT116 and HT29 
acts on the EGFR and not on downstream adaptors that are 
shared by many growth factor receptors. 

We next aimed to define the timescale on which the 
feedback operates. We treated HCT116 cells with the MER 
inhibitor for various times and then stimulated the pathway 
with EGF for 10 min and measured ART and ERR phosphor- 
ylation (Figure 5D}. Interestingly, already with a total 
inhibition time of 15 min, we observed an increase of ART 
phosphorylation, which increased further for longer pre- 
inhibition times, suggesting that the feedback operates as an 
integral feedback. Notably, after about 2h, also non-EGF- 
treated but MER-inhibited cells show a slight increase of ART 
levels, indicating an EGFR-independent cross talk operating 
on a longer timescale. Furthermore, the EGFR-dependent 
feedback shows that EGF treatment can partially restore the 
pre-inhibitor phospho-ERR level after ~4h of inhibition. 

EGFR inhibition prevents AKT activation by IVIEK 
inhibitors irrespective of mutation status 

We then used our model to predict combinatorial treatments 
that reduce ERR activity without activating ART. Our initial 
model was trained using only two inhibitors directly in the 
EGFR pathway, thus we decided to retrain the model for 
further inhibitors against EGFR (gefitinib) and RAF 
(sorafenib), see Figure 6A blue bars and Supplementary 
Figure S5. We then used the model to predict treatments of 
inhibitors of RAF or MER combined with EGFR or PI3R plus the 
combination of the latter two (red bars in Figure 6A} . For all 
these combinatorial perturbations, the model predicted a 
reduction in ART activity when compared with TGFa treatment 
only. Measurement of phospho-ART confirmed the predic- 
tions, except for combinations of MER/RAF inhibitors and 
PI3R inhibitors in HT29 cells (black bars in Figure 6A}. 

EGFR inhibition together with inhibition in IVIAPK 
signaling prevents growth of tumor cells 
irrespective of mutation status 

We next asked whether a combination therapy of MAPR 
inactivation together with an inhibitor against EGFR also 
stops proliferation of tumor cells. To assess this, we 
measured cell growth of HCT116 and HT29 cells with those 
four inhibitors alone, their vehicle control DSMO, and in 
selected combinations. Inhibition of the EGF receptor as well 
as inhibition of PI3R alone had no effect on proliferation in 
HT29 and HCT116 (Figure 6B}, and also combined applica- 
tion of both inhibitors did not alter proliferation either. This 
is in line with the notion that the oncogenes KRAS and BRAF 
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drive proliferation in these cells via MAPR signaling and do 
not require the EGFR to generate the pro-proliferative signal. 
We next investigated how manipulations downstream of 
their driver mutations alter proliferation in combination with 
PI3R inhibition. To inhibit the MAPR signaling pathway, we 
decided to use two different inhibitors in the two cell lines. 
Based on the differences in ERR-RAF feedback (Sturm et at, 
2010; Fritsche-Guenther et aU 2011), we decided to inhibit 
MAPR activity with the RAF inhibitor sorafenib in HCT116 
cells and the MER inhibitor AZD6244 in HT29 cells. This 
decision is further supported by the corresponding measure- 
ments of the double inhibition in Figure 6A, where those 
combinations showed the strongest effects. 

Application of sorafenib alone had no effect on growth of our 
RAS/PI3K mutant cell line model HCT116, suggesting that 
growth may depend on multiple redundant pathways. However, 
in combination with a PI3R inhibitor, the cells stop growing, 
indicating that MAPR and ART signaling redundantly control cell 
growth. As our model predicts that a combination of RAF and 
EGFR inhibitor prevents ART activation, we tested the EGFR 
inhibitor gefitinib together with sorafenib. In line with the 
model, this combinatorial treatment synergistically reduced 
growth as strong as the combination with PI3R inhibition. 

In Bi^AF mutant HT29 cells, MER inhibition blocked growth, 
and PI3R inhibition had no effect, confirming that HT29 cells 
depend solely on MAPR signaling for growth. The combina- 
tion of MER inhibitor and PI3R inhibitor led to a decrease of the 
cell index. Our model predicts that EGFR would also synergize 
with MER inhibition to block ART activation. In line with this 
prediction, EGFR inhibition had no effect when provided 
alone, but caused a decrease in cell index when cells were 
treated in combination with a MER inhibitor. 



EGFR synergizes with MEK inhibitors in vivo 

We sought to substantiate the effects of co-targeting of MER 
and EGFR in the DLD-1 colorectal xenograft model that 
harbors both KRAS^^^^ and PIKSCA^^"^^^ mutations. DLD-1 
tumors were estabhshed in nude mice and treated with 
erlotinib (dosed at a clinically relevant dose of 50mg/kg}, 
GDC-0973 (a potent MER-1/2 allosteric inhibitor, dosed both at 
1 and 5 mg/kg), or the combination of erlotinib + GDC-0973 at 
both dose levels. While erlotinib and the lower dose of GDC- 
0973 (1 mg/kg) were ineffective as single agents resulting in 1 
and 22% tumor growth inhibition (%TGI}, respectively, the 
combination demonstrated superior combination efficacy over 
either single agent alone with 36% TGI (Figure 6C, top panel). 
The 5 mg/kg dose of GDC-0973 demonstrated better single- 
agent activity with 42 % TGI, however again combination with 
erlotinib demonstrated superior tumor inhibition to either 
drug used alone with 60% TGI (Figure 6C, bottom panel). The 
combination of these drugs was well tolerated in mice with 
minimal weight loss (Supplementary Figure S6). 

Taken together, the results suggest that the EGF receptor is 
not required for maintaining the cellular phenotype in cells 
with RAS and i^AF mutations. However, once MAPR signaling 
is blocked, the EGFR increases ART activity and thus regains its 
key role in the network that was previously lost due to 
downstream mutations. 

© 2013 EMBO and Macmillan Publishers Limited 
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Figure 6 Model-derived combinatorial treatment options unveil effective combinations. (A) Model prediction of selected double inhibitions (red bars) and Luminex 
measurements (dark bars) for phosphorylation of AKT in HCT1 1 6 and HT29 cells in the presence of TGFa in log2 scale. Fits for single perturbations are shown on the left 
as blue bars. Error bars for modeled results represent s.d. of 100 noised simulations (see Supplementary Methods) (B) Cell proliferation (Xcelligence) of HCT1 16 and 
HT29 cells in response to single and double inhibition and control. Arrows mark time of treatment, and areas represent the maximal range of the growth curves, with dark 
lines reflecting the mean (n> = 3). For panels A and B, concentration of AZD6244 was 1 |iM. (C) Tumor growth of coloretal cancer cell DLD-1 cells transplanted into 
nude mice, which received a daily oral gavage of one of two different concentrations of MEK inhibitor GDC-0973, EGFR inhibitor eriotinib alone or combined. Error bars 
represent s.e.m. with n= 10. Source data for this figure is available on the online supplementary information page. 



Discussion 

A huge body of detailed mechanistic understanding about 
signaling processes has been accumulated w^ithin the last 
decades (Wang et al 2002; Oda et al 2005). Still it remains 
challenging to understand how^ a signaling netw^ork reacts 
v^hen targeted therapies are apphed because of strong cross 
talk betw^een pathw^ays and feedback loops (Kim et al 2007; 
Friday et al 2008; Sturm et al 2010; Fritsche-Guenther et al 
2011; Kholodenko et al 2012). In this study, w^e addressed this 
problem by quantifying MAPK/AKT signaling netw^orks in a 
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panel of colorectal cancer cell lines v^ith differing genetic 
background. We developed an experimental/computational 
pipeline that compiles mathematical models for individual cell 
lines based on systematic perturbation. Our modeling 
approach resides betw^een detailed mechanistic models 
(Schoeberl et al 2002) and more coarse-grain logical models 
(Saez-Rodriguez et al 2011, 2009} or even purely statistical 
descriptions (Tentner et al 2012} . The complexity of the model 
w^as chosen such that it could be parameterized v^ith the 
limited perturbation data set, but allow^ed to reveal netw^ork 
features such as feedback loops that simpler representations 
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cannot account for. MRA generally requires that the response 
of the system to perturbations can be modeled using linear 
equations, and the system is close to steady state. Thus, the 
modeling procedure can be used to interpret perturbation 
screens, but parameters should be interpreted in a phenom- 
enological rather than in a precise mechanistic way. Conse- 
quently, the procedure is helpful to interpret perturbation data 
and generate new hypotheses that can be subsequently tested, 
such as in a previous study of transient EGF/NGF signaling 
(Santos et aU 2007). In our study, the majority of parameters 
could be well estimated from the data. However, if there are 
strong uncertainties in the parameters, methods such as 
MCMC (Hastings, 1970} or the profile likelihood (Raue et al, 
2009} method can be readily apphed to model parameters or 
structural uncertainties. 

When we compared the signaling maps between the cell 
lines, we observed that the core signaling network was 
quantitatively very similar in all cells, although these cells 
were heterogeneous in their genetic constellation. This 
suggests that it is instrumental to build generic models of 
signaling despite genetic heterogeneity. Nevertheless, certain 
quantitative aspects differed strongly between cells, such as 
the strength of the response toward stimulation of the EGFR. 
We also found qualitative differences between cell lines, such 
as the loss of the ERK-RAF feedback in HT29 cells, which can 
be traced back to the BRAF ^^^^^ mutation (Friday et aU 2008} . 
Similar studies on larger cell line collectives may unveil further 
differences in network wiring due to the underlying mutations. 

Despite the diversity of mutations in the EGFR signaling 
network, we found a conserved strong feedback from ERK to 
EGFR in all five cell lines. At least in HCT116 cells, this 
feedback was detectable within ISmin but its strength 
increased further on timescales of hours, suggesting that this 
feedback operates as an integral feedback, which may have a 
role for signaling homeostasis (Yi et aU 2000; Prahallad et aU 
2012}. Different mechanisms of feedback regulation of the 
EGFR are known, which may all contribute. Transcriptional 
feedbacks such as MIG-6 (Yoon et aU 2012} or the Sprouty 
family (Mason etaU 2006; Halilovic etaU 2010} cannot account 
for the fast feedback regulation, but may be involved in later 
phases. Adaptors shared between receptors, such as SOS 
(Douville and Downward, 1997; Shankaran and Wiley, 2010} 
or Gabl (Yu et aU 2002} , are fast enough but most likely not the 
main feedback players, as the cross talk was only strong when 
EGFR was stimulated, but was not or only weakly present 
when HGF, FGF or IGF were apphed. Thus, the strongest target 
of this feedback is most likely the EGFR itself. Pancreatic 
cancer cells exhibited increased phosphorylation of the EGFR 
at Y1068, Y1045 and Y845 when MEK was inhibited (Gan et aU 
2010}. Phosphorylation of T669 by ERK affects EGFR turnover 
(Birtwistle et aU 2007}. Knockdown of CDC25A, known to be 
regulated by ERK (Wang et aU 2002}, did also lead to increased 
phosphorylation of EGFR at Y1068 (Prahallad et aU 2012}, 
suggesting that ERK changes CDC25A activity and by this 
regulates phosphorylation of EGFR (Wang et aU 2005}. 

A consequence of this feedback for targeted inhibition is 
that it leads to activation of AKT upon inhibition within 
MAPK signaling. Such feedback-mediated cross talk has 
been noted in many different tumors such as breast cancer 
(Mirzoeva et aZ, 2009; Lu et aU 2011}, prostate cancer (Gan 
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etaU 2010}, melanoma (Gopal etaU 2010}, gastric cancer (Yoon 
et aU 2009} and in colorectal cancer (Prahallad et aU 2012}. 
Increased phosphorylation of AKT may cause drug resistance, 
as AKT activity stimulates survival and migration. Further- 
more, AKT shares a complex network of transcription factors 
with ERK (Stelniec-Klotz et aZ, 2012} and both paths converge 
on key proteins important for cellular function such as cyclin- 
Dl for growth control (Halilovic et aU 2010} and Bad for 
apoptotic regulation (She et aU 2005}, explaining the need to 
switch off both pathways. 

Our model allowed devising combinatorial therapies that 
block ERK activation and at the same time prevent rise in 
AKT activity. While PI3K inhibitors may be used to block AKT 
signaling, EGFR inhibition can also prevent strong AKT 
activation, irrespective of whether RAF, RAS or PI3K was 
mutated. In line with this prediction, we found approximately 
the same synergistic reduction of growth for PI3K/MAPK 
combination than for EGFR/MAPK inhibition in vitro in two 
cell line models. Therefore, while upregulation of AKT by 
MAPK inhibition can be successfully blocked by PI3K or mTOR 
inhibition (Balmanno et aZ, 2009; Mirzoeva et aZ, 2009; 
Aksamitiene etaU 2010} in colorectal cancer models, upstream 
inhibition of EGFR may be similarly potent with possibly less 
side effects. Using a RAS/PI3K mutant colon cancer xenograft 
model, we confirmed that combinatorial therapy with inhibi- 
tors against MEK and EGFR is also an efficient therapy in vivo. 
This extends previous findings showing that BRAF inhibitors 
synergize with EGFR inhibition in a colorectal cancer 
xenograft model with BRAF mutation (Prahallad et aZ, 2012}. 
For i^AS-mutated tumors, so far no targeted therapy is 
available (Baines et aZ, 2011; Ward et aZ, 2012}, and a mutation 
in RAS precludes EGFR-directed interventions. Our results 
suggest that i^AS-mutated tumor cells can be successfully 
treated by EGFR inhibitors if provided together with MEK or 
RAF inhibitors. 

While our model can predict successful combinatorial 
treatments, it has limitations. It cannot account for combina- 
tions that require sequential appHcation of drugs (Lee et aZ, 
2012}, and fails to capture resistance due to tumor-stroma 
interactions (Sebens and Schafer, 2012} . It is likely that in other 
cell types different combinations of drugs may be more 
successful, as the role of specific feedbacks can be different in 
different cell types, and can even switch between positive and 
negative effects depending on receptor expression (Birtwistle 
et aZ, 2007}. 

Tumor evolution is one of the major causes for eventual 
relapse (Iwasa et aZ, 2006}. For example, relapse after long- 
term treatment with the anti-EGFR agents Panitumumab in 
KRAS wild-type colorectal carcinoma (Amado et aZ, 2008; 
Karapetis et aZ, 2008} is often caused by selection for mutations 
downstream of EGFR (Diaz et al, 2012; Misale et aZ, 2012} . The 
combinatorial treatment predicted by our model may thus be 
advisable even for EGFR-addicted tumors, as it will counteract 
selection for additional mutations in RAS or RAF. 

In colon cells, the EGFR receptor is very potent, but 
depending on tissue and mutational status, other receptors 
may be more important in stimulating the network. For 
example, c-Met, that is feedback regulated and signals to both 
ERK and AKT, was shown to mediate resistance to BRAF 
inhibition in melanoma (Wilson et aZ, 2012}, suggesting c-Met 
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as the most effective co-target for this particular situation. 
Hence, our results further highlight that downstream muta- 
tions such as in RAS and BRAF do not necessarily invalidate 
upstream drug treatment if provided in combination with a 
suitable downstream inhibitor. 

Materials and methods 
Model construction and evaluation 

The general modeling approach is depicted in Figure 3A. First, we 
derived a starting network from the literature. We then developed a 
modeling framework that extends our previous algorithm (Stelniec- 
Klotz et al, 2012) to estimate a quantitative map of the signaling 
network using an approach derived from MRA (Kholodenko, 2000; 
Bruggeman et al, 2002; Friday et al, 2008; Cirit et al, 2010; Sturm et al, 
2010; Fritsche-Guenther et al, 2011). Briefly, MRA links the log-fold 
change after perturbation (called global response coefficient R] to a 
matrix corresponding to link strengths (local response matrix r) by 
inversion. We adjust the equation linking R and r for two qualitatively 
different perturbations (stimulation and inhibition) . Stimulations are 
modeled as entries in the perturbation vector and inhibitions impair 
the signal flow of these stimulations. In addition, we allowed inhibitors 
to affect the base level of the signaling of the inhibited node, which is 
incorporated in the perturbation vector as stimulus. As the experi- 
mental design allowed measuring and perturbing only selected nodes, 
many parameters were not identifiable. We developed an algorithm 
using Gaussian elimination that reparameterized the model such that 
all parameters became identifiable. These new parameters were then 
determined by maximizing the likelihood using a Levenberg- 
Marquardt optimization algorithm using the best fit of 10 000 random 
initializations of parameters. We verified that this initial parameter 
scan results in unique parameters, as running the procedure 100 times 
on the same data, each time initializing the random number generator 
with a different number resulted in identical parameter sets. This 
network was then pruned by using a greedy hill climbing approach in 
which we removed edges that after model refitting did not significantly 
decrease the goodness of fit (Likelihood ratio test, with significance 
threshold of 0.05). The implementation of this algorithm was 
essentially as described earlier (Stelniec-Klotz etal, 2012) but extended 
by the non-identifiability analysis and combinatorial perturbations. To 
extend the model for further inhibitors, we performed maximum 
likelihood estimates for the novel parameters, while keeping the others 
constant. Error bars for predictions were generated by 100 times fitting 
the model to data with added noise according to the error model. The 
program was written in the programming language C^, and 
symbolic matrix operations were conducted with the GiNaC library 
version 1.6 (http://www.ginac.de/). An extensive description of the 
algorithm and results of a stepwise reduction can be found in 
Supplementary Method 1 and Supplementary Figure S2, respectively. 



Cells and cell culture 

Human colorectal cancer cell lines SW480, SW403, HCT116, RKO and 
HT29 were obtained from the ATCC (American Type Culture 
Collection, UK). LIM1215 were kindly provided by Professor John 
Mariadason (Ludwig Institute for Cancer Research Austin Hospital, 
Melbourne), Australia. All cell lines were maintained in DMEM 
(Dulbecco's modified Eagle's medium, Lonza) supplemented with 
10% fetal calf serum, 1% ultraglutamine and 1% penicillin/ 
streptomycin and incubated in a humidified atmosphere of 5 % C02 
in air at 37 °C. 



Reagents 

The following inhibitors where used in the various assays: MEK 
inhibitors AZD6244 (0.1 |iM unless otherwise specified, Selleck 
Chemicals LLC) and GDC-0973 (Genentech), RAF inhibitor sorafenib 
(10|iM, LC Laboratories), PI3K inhibitor LY294002 (10 |iM, Axxora 
Deutschland), EGFR inhibitors gefitinib (10|iM, Cayman Chemicals) 



and erlotinib (Genentech), GSK3B inhibitor SB216763 (5 |iM, Sigma 
Aldrich), IKK-a/P inhibitor BMS345541 (25 |iM, Sigma Aldrich), and 
IKK-B inhibitors PS1145 (10|iM, AxonMedChem) and PHA408 
(100 nM, AxonMedChem). The solvent control was DMSO (equal |iM 
to each inhibitor) . 

We used the following ligands (all Peprotech): TGF-a (0.01 |ig/ml), 
IGF-1 (0.1|ig/ml), EGF (0.025 jig/ml), FGF2 (0.005 |ig/ml) and HGF 
(0.05 |ig/ml) with 0.01 % BSA in PBS as solvent. 



Luminex assays 

After treatment and incubation, lysates were collected according to 
supplier's protocol and analyzed with the Bio-Plex Protein Array 
system (Bio-Rad, Hercules, CA) using beads specific for P-AKT (S473), 
P-ERKl/2 (Thr202/Tyr204/Thrl85/Tyrl87), P-ERK2 (Thrl85/ 
Tyrl87), P-GSK3a/B (S21/S6), P-IGF-IR (Tyrll31), P-IRSl (S636/ 
S639), P-IKBa (S32/S36), P-MEKl (S217/S221), and P-p70S6K 
(Thr421/S424) according to the manufacturer's instructions. Briefly, 
samples were washed with PBS and lysed with cell lysis buffer 
(Bio-Rad). Lysate protein concentration was determined with BCA 
(bicinchoninic acid) method. The beads and detection antibodies were 
diluted 1:5 or 1:3. For data acquisition, the Bio-Plex Manager software 
was used. 



RNA isolation and quantitative RT-PCR analysis 

RNA was isolated from HCT116 cells treated with BMS345541 using the 
RNeasy-mini-kit (Qiagen) according to the supplier's protocol. To 
obtain cDNA from RNA, the high-capacity cDNA reverse transcription 
kit (Applied Biosystems) was used. Synthesis of double-stranded DNA 
during the PGR cycles was visualized with TaqMan gene expression 
assays FAM-dye labeled (gene of interest: EGRl (Hs_00152928_ml), 
cFos (Hs_00170630_ml)) or VIC-dye labeled for loading control PGKl 
(Hs_943178_gl) and TaqMan gene expression master mix (Applied 
Biosystems). Quantitative real-time PCR (qrt-PCR) analysis was 
performed using a StepOnePlus 96-well format Light-Cycler apparatus 
(Applied Biosystems). Experiments were run and analyzed with the 
StepOne 2.0 software. The data were analyzed quantitatively by 
measuring the threshold cycles (Ct) . Cj values where normalized first 
by the Cyof the internal control PGKl ( - ACt) and second by the ACt 
value of the zero time point ( - AA Ct), which are interpreted as log2 
fold changes. 



Immunoblotting 

Protein extracts of cells were prepared as described for Bioplex 
analysis. Blotting procedure and materials were as previously 
described (Fritsche-Guenther et al, 2011; Stelniec-Klotz et al, 2012). 
The following primary antibodies were used: rabbit anti-human 
P-p70S6 (Thr389, Cell Signaling Technology, 1:500) or mouse anti- 
human GAPDH (Ambion, 1:12 500). Membranes were scanned using 
Li-COR Odyssey. The signals were quantified using Odyssey software. 



Tumor and body weight measurements 

The DLD-1 colorectal tumor cell line was inoculated subcutaneously 
on the flank of female athymic nude (nu/nu) mice (Harlan 
Laboratories) and tumors were allowed to establish to a size of 
125-250 mm^. Animals were grouped out into six treatment groups 
[n = 10 per group) based upon tumor volume and treatment was 
initiated with vehicles (MCT (methylcellulose + 0.2% Tween 80) and 
15% HBCD (15% hydroxypropyl-beta-cyclodextrin)), erlotinib 
(50mg/kg, daily oral gavage formulated in 15% hydroxypropyl-beta- 
cyclo dextrin), GDC-0973 (1 or 5mg/kg, daily oral gavage formulated 
in MCT), or the combination of erlotinib + GDC-0973 at both doses. 
Tumor volumes were determined using digital calipers (Fred V. Fowler 
Company, Inc.) using the formula (LxWxW)/2. %TGI was 
calculated as the percentage of the area under the fitted curve (AUG) 
for the respective dose group per day in relation to the vehicle, such 
that % TGI = 100 X 1 — (AUCtreatment/day) / (AUCvehicie/day) . Curve 
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fitting was applied to log2-transformed individual tumour volume data 
using a linear mixed-effects model using the R package nlme, version 
3.1-97 in R v2. 12.0. Body weights and gross observations on animal 
welfare were also recorded throughout the study. All experimental 
procedures conformed to the guiding principles of the American 
Physiology Society and were approved by Genentech's Institutional 
Animal Care and Use Committee. 



Proliferation assay 

For proliferation studies in real time, the XCelligence RTCA SP 
instrument was used. HCT116 and HT29 cells were plated 24 h before 
treatment with inhibitors or controls. Over the time of measurement, 
the system records the cell index (ci), which reflects the cell 
attachment to the electrodes and is proportional to the number of 
cells. Each treatment was measured at least in triplicates. As the initial 
cell number is variable and growth of the cell population is 
exponential, averaging the growth curves is not reasonable. Instead, 
we calculated the logarithmised and normalized cell index log2(ci(0/ 
ci[t of inhibitor application}} and took the average of these values. 



Sequencing and detection of single-nucleotide 
variations (SNVs) 

We performed targeted sequencing by the Ion AmpliSeq Cancer Panel 
1.0 and the Ion PGM sequencer covering 13.7-kb-coding sequence 
from 46 cancer-relevant genes (Supplementary Table ST2} and 739 
COSMIC-annotated SNP positions. Genomic DNA was extracted from 
pelleted CRC cell lines. The extracted DNA was quantified and the 
quality estimated by the Qubit 2.0 Fluorometer (Qubit dsDNA HS 
Assay Kit, Life Technologies} and Bioanalyzer 2100 (High Sensitivity 
DNA Kit, Agilent} . The amplicon library was prepared according to the 
manufacturer's protocols using 12-30 ng of DNA and loaded on an Ion 
314 Chip (Life Technologies} to yield at least 10 Mb. To increase the 
efficiency, we spun the chip for 2 min before and after turning around 
the chip by 180° for each loading. lonTorrent Suite, Version 2.1 Variant 
Caller (Life Technologies} was used for variant calling, with positions 
with minor allele frequencies below 10 % were defined homozygotic. 
SNVs were filtered as described in supplement. 



Supplementary information 

Supplementary information is available at the Molecular Systems 
Biology website (www.nature.com/msb}. 
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